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Abstract. In this work, we present a new approach to the construction of variational 
integrators. In the general case, the estimation of the action integral in a time interval 
[qk,qk+i] is used to construct a symplectic map (qk,qk+i) — * (<7fc+i> Qk+s)- The basic 
idea here, is that only the partial derivatives of the estimation of the action integral 
of the Lagrangian are needed in the general theory. The analytic calculation of these 
derivatives, give raise to a new integral which depends not on the Lagrangian but on 
the Euler-Lagrange vector, which in the continuous and exact case vanishes. Since this 
new integral can only be computed through a numerical method based on some internal 
grid points, we can locally fit the exact curve by demanding the Euler-Lagrange vector 
to vanish at these grid points. Thus the integral vanishes, and the process dramatically 
simplifies the calculation of high order approximations. The new technique is tested 
for high order solutions in the two-body problem with high eccentricity (up to 0.99) 
and in the outer solar system. 
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1. Introduction 

It is well known that the dynamics of seemingly unrelated conservative systems in 
mechanics, physics, biology, and chemistry fit the Hamiltonian formalism ([]]). Included 
among these are particle, rigid body, ideal fluid, solid, and plasma dynamics. An 
important property of the Hamiltonian flow or solution to a Hamiltonian system is 
that it preserves the Hamiltonian and the symplectic form (see, for example, [2]). A 
key consequence of symplecticity is that the Hamiltonian flow is phase-space volume 
preserving (Liouvilles theorem). Since analytic expressions for the Hamiltonian flow are 
rarely available, approximations based on discretization of time are used. A numerical 
integration method which approximates a Hamiltonian flow is called symplectic if it 
discretely preserves a symplectic 2-form to within numerical round off ([3j HJ E]) and 
standard otherwise. By ignoring the Hamiltonian structure, a standard method often 
introduces spurious dynamics. This effect is excellent illustrated in [6] where a spherical 
pendulum with symmetric perturbation is integrated using implicit Euler and a 
symplectic method (variational Euler). The computation of a Poincare section shows 
that the standard method, despite the order of magnitude in integration time-step 
size, artificially corrupts phase space structures by exhibiting a systematic drift in the 
invariant tori, whereas variational Euler preserves them. Moreover, in systems that are 
non-integrable, symplectic integrators often perform much better even compared with 
projection methods, where the calculated values are projected in the manifold defined by 
the constant symplectic structure. Finally, as it shown in [7], in many body problems, 
the symplectic integrators perform increasingly better than standard methods as the 
number of bodies increases. 

Symplectic integrators can be derived by a variety of ways including Hamilton- 
Jacobi theory, symplectic splitting, and variational integration techniques. Early 
investigators, guided by Hamilton- Jacobi theory, constructed symplectic integrators 
from generating functions which approximately solve the Hamilton-Jacobi equation 
(P HJ [5]). The symplectic splitting technique is based on the property that symplectic 
integrators form a group, and thus, the composition of symplectic-preserving maps is 
also symplectic. The idea is to split the Hamiltonian into terms whose flow can be 
explicitly solved and then compose these individual flows in such a fashion that the 
composite flow is consistent and convergent with the Hamiltonian flow being simulated 
as it is explained in detail in [8]. On the other hand, variational integration techniques 
determine integrators from a discrete Lagrangian and associated discrete variational 
principle. The discrete Lagrangian can be designed to inherit the symmetry associated 
with the action of a Lie group, and hence by a discrete Noethers theorem, these methods 
can also preserve momentum invariants (for a discussion of the above statements, see 

[ID- 

Variational integration theory derives integrators for mechanical systems from 
discrete variational principles (0, [101 HH H2J EE]). Variational principles have been 
succesfully applied to partial differential equations and to stochastic systems as well 
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see, for example, [14J). In the general theory, discrete analogs of the Lagrangian, 
Noethers theorem, the Euler-Lagrange equations, and the Legendre transform can be 
easily obtained [T5l [TBI [T7j. Moreover, variational integrators can readily incorporate 
holonomic constraints (via Lagrange multipliers) and non-conservative effects (via their 
virtual work) ( [TT1 fl~2]). The algorithms derived from this discrete principle have been 
successfully tested in infinite and finite-dimensional conservative, dissipative, smooth 
and non-smooth mechanical systems (see [18] [191 f20j and references therein). Recently, 
variational principles have been applied to particle mesh methods [21] and to fractional 
stochastic optimal control [22]. In the general approach, as it will be presented in 
section [21 the variational principle to apply is the estimation of the action integral in a 
time interval [ioj^i] as a smooth function of the edges of the interval <7o>(Zi- Since any 
sufficiently smooth and non-degenerate function S(qo, qi) generates via 

dS(q ,qx) dS{q Q ,q 1 ) 

P° = q >Pi = o 

dq dq 1 

a symplectic map (ftbPo) ~~ ¥ (<?i>Pi) ([23]), the estimated action integral can be used to 
develop a discrete analog to the Euler-Lagrange equations. 

The accuracy is not the terminus of the application of variational integrators, but 
rather their ability to discretely preserve essential structure of the continuous system 
and in computing statistical properties of larger groups of orbits, such as in computing 
Poincare sections or the temperature of a system (see, for example, [20], [24]). On the 
other hand, high accuracy can be obtained using special designed methods as it is 
explained in [25] [26]. These methods, although they produce high order estimations 
of the positions and the momenta, computationally become very heavy as the order 
increases. This is due to the large number of parameters that have to be calculated 
in each step as a solution to a non-linear system, depending on the nature of the 
Lagrangian. The total set of equations, consists of partial derivatives of the action 
integral and a set of variational equations that determine the internal points, necessary 
for high order methods. 

In this work, we present a new approach to the construction of variational 
integrators. The basic idea, is that only the partial derivatives of the estimation of 
the action integral of the Lagrangian are needed in the general theory. The analytic 
calculation of these derivatives, give raise to a new integral which depends not on the 
Lagrangian but on the Euler-Lagrange vector, which in the continuous and exact case 
vanishes. Since this new integral can only be computed through a numerical method 
based on some internal grid points, we can locally fit the exact curve by demanding the 
Euler-Lagrange vector to vanish at these grid points. Thus the integral vanishes, and 
the process dramatically simplifies the calculation of high order approximations. 

2. Discrete Variational Mechanics 

The well known least action principle of the continuous Lagrange - Hamilton Dynamics 
can be used as a guiding principle to derive discrete integrators. Following the steps of 
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the derivation of Euler-Lagrange equations in the continuous time Lagrangian dynamics, 
one can derive the discrete time Euler-Lagrange equations. For this purpose, one 
considers positions go and gi and a time step h G R, in order to replace the parameters 
of position q and velocity q in the continuous time Lagrangian L(q,q,t). Then, by 
considering the variable h as a very small (positive) number, the positions go an d qi 
could be thought of as being two points on a curve (trajectory of the mechanical system) 
at time h apart. Under these assumptions, the following approximations hold: 

qo ~ q(0) , qi « q(h) , 

and a function L d (q , g l5 h) could be defined known as a discrete Lagrangian function. 

Many authors assume such functions to approximate the action integral along the 
curve segment between qo and qi, i.e. 

L d (q ,q 1 ,h)= / L(q,q,t)alt (1) 
Jo 

Furthermore, one may consider the very simple approximation for this integral given 
on the basis of the rectangle rule described in [TT] . According to this rule, the integral 
J Q T Ldt could be approximated by the product of the time-interval h times the value of 
the integrand L obtained with the velocity q replaced by the approximation (q 1 —q )/h: 
The next step is to consider a discrete curve defined by the set of points {q k } k= o, and 
calculate the discrete action along this sequence by summing the discrete Lagrangian of 
the form L d (q k , qk+i, h) defined for each adjacent pair of points (q^, q k +\)- 

Following the case of the continuous dynamics, we compute variations of this action 
sum with the boundary points go an d g^v held fixed. Briefly, discretization of the action 
functional leads to the concept of an action sum 

n-l 

Sd{ld) = ^2,L d (q k -i,q k ), J d = (go, •••,9n-l) G Q n (2) 

k=l 

where L d : Q x Q — > R is an approximation of L called the discrete Lagrangian. 
Hence, in the discrete setting the correspondence to the velocity phase space TQ is 
QxQ. An intuitive motivation for this is that two points close to each other correspond 
approximately to the same information as one point and a velocity vector. The discrete 
Hamilton's principle states that if 7^ is a motion of the discrete mechanical system then 
it extremizes the action sum, i. e., 5S d = 0. By differentiation and rearranging of the 
terms and having in mind that both g and q^ are fixed, the discrete Euler-Lagrange 
(DEL) equation is obtained: 

D 2 L d (qk-i, Qk, h) + DiL d (q k , q k +i, h) = (3) 

where the notation DiL d indicates the slot derivative with respect to the argument of 
L d . 

We can define now the map $: Q x Q Q x Q, where Q is the space of generalized 
positions q, by which 

D x L d o $ + D 2 L d = (4) 
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which means that $(q k -i,q k ) — (<?fc?<7fc+i)- Then, if for each q e Q, the map 
D 1 L c [(q, q) : T q Q — > T*Q is invertible, then DiL d : Q x Q — ► T*Q is locally invertible 
and so the discrete flow defined by the map $ is well defined for small enough time steps 
(see [27] for details). Moreover, if we define the fiber derivative 

FL d :QxQ^ T*Q (5) 

and the two-form u on Q x Q by pulling back the canonical two-form VLqan = dq % A dpi 
from T*Q to Q x Q: 

u = FL* d (Q CAN ) (6) 

The coordinate expression for u is 

dq{dq J k+1 

and can be easily proved that the map $ preserves the symplectic form uj (two different 
proofs are presented in [28] and [12]). Finally, assuming that the discrete Lagrangian 
is invariant under the action of a Lie group G on Q and £ G g, the Lie algebra of 
G, by analogy with the continuous case, we can define the discrete momentum map 
Jd ■ Q x Q -> by 

(</d(<7fc, gfe+i), := (D a L d (q k , q h+1 , € Q (q k )} (8) 

It can be proved that the map <!> preserves the momentum map J d [12J. 

In a position-momentum form the discrete Euler-Lagrange equations ([3]) can be 
defined by the equations below 

p k = - D x L d (q k , q k +i, h) 

p k+1 = D 2 L d (q k , q k+ x, h) (9) 
3. Local Path Fitting 

The system in eqj9]can be considered now as a numerical one-step method {q k ,Pk) 
(q k+ i,p k+ x). The level of the accuracy in estimation of the integral 

L d {q k ,q k +x,h) ~ / L(q,q,t)dt (10) 
Jt k 

fully characterizes the accuracy of the method. But as we can see from eqJHJ only the 
derivatives of L d are needed. Thus, consider a parameter A. Then 

dL d (q k ,q k+ i,h ) _ Tf _ u N . u ^ u ^dh , f tk+1 dL(q,q,t) 
d\ 

But, 



x ,*dh f tk+1 dL(q,q,t) , 
L{q{t k+1 ),q{t k+1 ),h)— + J^ } dt (11) 



9L(q,q,t) _ dL(q,q,t) dq dL(q,q,t) dq 

d\ dq d\ dq dX 1 ' 

or, changing the order of derivation in the second term of the right hand part, 

dL(q,q,t) = dL(q,q,t) dq dL(q,q,t)d§^ 

dX dq dX dq dt 1 } 
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Integrating by parts now, we get 
dL(q,q,t) 



dh 



dX 

tk+l 



L(q(t k+1 ),q(t k+1 ),h) — + 
( dL(q,q,t) al ( dL(q,q,t) 



dq_ df dL(q, q,t) dq *fc+i 
d\ dq dX 



t \ dq dt \ dq 

Now, instead of estimating the integral in eqJTUJ we only have to estimate 



f dL(q,q,t) _ d_ f dL(q,q,t) \\ dq 
V dq dt\ dq J J dX 



(14) 



(15) 



We can use any quadrature rule that is based on a set of S + 1 grid points at times 
{t k + c^h.j = 0, 1, ... , S}. On the other hand, if we demand that the Euler-Lagrange 
equation 

dL(q,q,t) _d_( dL(q,q,t) \ = Q 
dq dt \ dq J 

holds at these grid points, then J = and 
dL(q,q,t) 



(16) 



dX 



Tl ( . x v . \ 7 \ dh dL(q,q,t) dq 



tk 



;i7) 



The set of equations eq|9] are now given: 



Pk 



Pk+i 



9L(q,q,t) dq 



c fe+i 



dq dq k 
dL(q,q,t) dq 



dq dq k+1 

The above set of equations is consistent with the variational principles. Consider the 
S + 1 grid points at times P = t k + c^h, with c° = and c s = 1 the edge points q k , q k +i- 
For the internal points we have 
dL(q, q, t) dq 
dq dqi t k 

where g J are the internal points, because the curve is fixed at its endpoints and thus 



(19) 



dLd 

dqi 







(20) 



Let us now work an exact example. Consider the Lagrangian of the harmonic 
oscillator with unity frequency 



L = -q 2 



—Q 



and let q{t) in the interval [t k ,t k+ i] is given 



q(t) = q k [l 



t 



t 



1 



+ y 



l 



(21) 



(22) 



where h = t k+ \ — t k and x, y are free parameters. We demand now that the Euler- 
Lagrange equation holds at points t = t k and t = t k+1 . Any quadrature now for the 
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Time(x 2n) 



9 10 



Figure 1. The exact orbit (solid line) for the harmonic oscillator with unity frequency 
and the calculated points (o) for the first 10 periods using an orbit satisfying the Euler- 
Lagrange equation at the edge points. 



calculation of the action integral based only on edge points will give the set of equations 
eqflBl The parameters x, y are easily calculated 

h 2 ( , q k +i\ 
h 2 

y = -r (Qk+i - ?fc) 



Then, the method eqJT8l gives 

6hq k + q k (6- 2h 2 ) 



Qk+l 



Qk+l 



Q + h 2 
Qk+i — Qk h 



h 



, , Qk+l 

s{ qk + — 



in fig. [U the calculated positions are plotted for the first 10 periods using a step size 
h = 0.01. 



4. Local Path fitting using Bernstein basis polynomials 

The n + 1 Bernstein basis polynomials of degree n are defined as 



71 



J 



xi(l-x) n ~i , j = 0,1, 



71 



(23) 



where the binomial coefficient 
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The Bernstein basis polynomials of degree n form a basis for the vector space II n 
of polynomials of degree n. Some properties of the Bernstein basis polynomials are 
summarized in the following: 

• b j:Tl (x) = if j < or j > n. 

• bj tTl (0) = 5j t o and bj >n (l) = 5j jn , where 5 is the Kronecker function. 

• 6j,„(l -x) = b n -j, n (x) 

• b 'j,ni X ) = 71 ( 6 J-l,n-l(») - hn-l(x)) 

• The Bernstein basis polynomials of degree n form a partition of unity, i.e. 

E"= W^) = 1 

• A Bernstein polynomial can always be written as a linear combination of 
polynomials of higher degree, i.e. 

bj, n -i(x) = ^b j}n {x) + 3 -^b j+hn {x). 

The most important property of Bernstein basis is their capability to approximate 
continuous functions. Let f(x) be a continuous function on the interval [0, 1]. If we 
consider now the Bernstein polynomial 



Bn(f)(x) 



3=0 



h 3,n{x) 



then 



lim B n {f){x) = f{x) 



(25) 



(26) 



uniformly on the interval [0,1]. To prove this, suppose if is a random variable 
distributed as the number of successes in n independent Bernoulli trials with probability 
x of success on each trial; in other words, K has a binomial distribution with parameters 
n and x. Then we have the expected value E(K/n) = x. Applying the weak law of 
large numbers of probability theory 

K - - (27, 



lim P 

n-+oo 



— X 



n 



> 5 







for every 5 > 0. Because /, being continuous on a closed bounded interval, must be 
uniformly continuous on that interval, we can infer a statement of the form 



lim P 

n— >oo 



Consequently 



lim P 

n—*oo 

lim P 



K 

n 

K 

n 

~) 

n I 



> e 



K 

n 

f'l 

n 



0. 



+ 



(28) 



Elf 



fix) 



V 



> e/2 



E^f 

> e/2^j + 
= 



n 



-fix) 



> e 



Local Path Fitting: A New Approach to Variational Integrators 



9 



And so the second probability above approaches as n grows. But the second probability 
is either or 1, since the only thing that is random is K, and that appears within 
the scope of the expectation operator E. Finally, observe that E(f(K/n)) is just the 
Bernstein polynomial B n (f)(x). 

A more general statement for a function with continuous A;-th derivative is 



B n (fY k) \\oo < ^||/ (fc) |U and ||/« - B n (f) 



(*)| 







(29) 



where additionally = (l — (l — ^) • • - (l — ^) is an eigenvalue of B n and the 
corresponding eigenfunction is a polynomial of degree k. 

Consider now the Lagrangian L(q, q, t) of a given system and the state vector (qk, Pk) 
at a given time t k . Let x\ j = 0, 1, . . . , S a set of S + 1 parameters and the polynomial 

' ' ' " fh S * (30) 



q(t) = J2^hs( 

3=0 ^ 



h 



which estimates the position at the interval [tk, tk+i] with h = t k+ i—tk- Since q(tk) = x° 
and q(tk+i) = x s , we have 



q(t) = q k b ,s [ J + ^l^b^s 

' 1=1 



h 



t — t k \ , ( t — t k 
+ Qk+ibs,s 



and 



5-1 



m = E 



S ( ft-t k 

ft+l^ I 65-1,5-1 I — 



t — t k 
h 



bs-i, 



h 

t-t k 
h 



h 



(31) 



t — t k 
h 



+ 



(32) 



Replacing now the position and velocities in the Lagrangian and demanding the Euler- 
Lagrange equation to hold in grid points, we get 



dL 

Pk = -jrr 
oq 



dL 
dq 

Pk+i 



t=t k 

t=t k +cih dt \ dq 
dL 



t=t k +cih 



0, j = l,2,...,S-l 



dq 



t=tfc+i 



(33) 



The above system is solved for (x j , p k +i), , j = 1,2, 
q k+ i = X s . 



, S, and gives the next point 



5. Numerical Tests 

5.1. The 2-body problem 

We now turn to the study of two objects interacting through a central force. The most 
famous example of this type, is the Kepler problem (also called the two-body problem) 
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Figure 2. The relative error in energy during one period for eccentricity e = 0.5, step 
size h — 0.05 and for different values of 5 (number of intermediate points). (I) for 
5 = 3, (0) for 5 = 5, (x) for 5 = 7, (□) for 5 = 9 and (o) 5 = 11. 



that describes the motion of two bodies which attract each other. In the solar system 
the gravitational interaction between two bodies leads to the elliptic orbits of planets 
and the hyperbolic orbits of comets. 

If we choose one of the bodies as the center of our coordinate system, the motion 
will stay in a plane. Denoting the position of the second body by q = (qi,q2) T , the 
Lagrangian of the system takes the form (assuming masses and gravitational constant 
equal to 1) 



£(q,q>*) = ^q^q + -- 
2 q 



The initial conditions are taken 



q 



q 



o. 



1 - e 



(34) 



(35) 



where e is the eccentricity of the orbit. In the first experiment, we take the eccentricity 
e = 0.5 and step size h = 0.05 and plot the relative error in the energy during one period 
for several number of intermediate points S. The results are shown in fig. [2] where it 
is clear that the order of the approximation is increased with increasing number of 
intermediate points. 

In the second experiment, we take the eccentricity e = 0.99 and use an adaptive 
time step control in order to keep the relative error in energy smaller than 10~ 7 . Tabled] 
shows the number of integration steps needed for one period. The order of approximation 
is again clearly increases with increasing S as the mean step size needed for the same 
error in energy increases from 1.8 • 10~ 3 for S = 5 to 0.1 for S = 12. Finally, in the last 
experiment we integrate the two-body problem with eccentricity e = 0.99 for 10 4 periods 
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Table 1. Number of integration steps 



s 


No of Steps 


3 


> 10 4 


4 


> 10 4 


5 


3526 


6 


460 


7 


421 


8 


181 


9 


142 


10 


112 


11 


98 


12 


59 





Figure 3. Long term integration of the 2-body problem with eccentricity e = 0.99 
for 10 3 periods with S = 12 intermediate points and a step size which is adaptively 
calculated in order to keep the relative error in energy less than 1CP 7 . In the first plot, 
the calculated positions (+) and the exact ones (o) are plotted for the last period. In 
the second plot, the relative error in angular momentum is plotted. 



in order to check the long term behavior of the method. The number of intermediate 
points is S = 12 and the step size is adaptively controlled in order to keep the relative 
error in energy less than 10 -7 . Fig. [3] shows the results. In the upper sub-figure, we 
plot the position during the last period along with the exact solution while in the other, 
the relative error in angular momentum. 
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Time (in days) Time (in days) 

Figure 4. Integration of the outer solar system for 10 6 days with step size h = 50 
days and S — 6 intermediate points. The planets (J-upiter, S-aturn, U-ranus, N-eptune 
and P-luto) follow constant orbits, while the relative error in energy is ~ 10~ 7 , the 
maximum relative error in momentum is less than 10 -10 and finally the relative error 
in angular momentum is ~ 10 . 

5.2. The 5- Outer Planet System 

The next problem concerns the motion of the five outer planets relative to the sun. 
The problem falls in the category of the N-Body problem which is the problem that 
concerns the movement of N bodies under Newton's law of gravity. The Lagrangian of 
this system is 

Ltf, <f , t) = \ £ m < iff <t + G ± £ |S=^T < 36 > 

3=0 j=l k=0 11 ^ ^ 11 

where G is the gravitational constant, m? is the mass of body j and q 1 , q 1 are the vectors 
of the position and velocity of body %. In [23] the data for the five outer planet problem is 
given (these data are summarized in table [2]). Masses are relative to the sun, so that the 
sun has mass 1. In the computations the sun with the four inner planets are considered 
one body, so the mass is larger than one. Distances are in astronomical units, time is in 
earth days and the gravitational constant is G = 2.95912208286 • 10~ 4 . The Lagrangian 
fl36l) has been integrated for t e [0, 10 6 ] with step size h = 50 days and using S = 6 
intermediate points. The results are shown in fig. HI The planet orbits are stable, the 
maximum relative error in energy is ~ 10 -7 , the maximum relative error in momentum 
is less than 10~ 10 and finally the relative error in angular momentum is ~ 10~ 9 . Note 
here, that errors in angular momentum are caused by the round off error produced by 
the solution of the non-linear system used to calculate the intermediate points. 
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Table 2. Initial Data for the 5-outcr planet problem 



Planet 


Mass 


Initial Position 


Initial Velocity 


Sun 


1.00000597682 




























Jupiter 


0.000954786104043 


-3.5023653 


0.00565429 






-3.8169847 


-0.00412490 






-1.5507963 


-0.00190589 


Saturn 


0.000285583733151 


9.0755314 


0.00168318 






-3.0458353 


0.00483525 






-1.6483708 


0.00192462 


Uranus 


0.0000437273164546 


8.3101420 


0.00354178 






-16.2901086 


0.00137102 






-7.2521278 


0.00055029 


Neptune 


0.0000517759138449 


11.4707666 


0.00288930 






-25.7294829 


0.00114527 






-10.8169456 


0.00039677 


Pluto 


1/(1.3- 10 s ) 


-15.5387357 


0.00276725 






-25.2225594 
-3.1902382 


-0.00170702 
-0.00136504 



6. Conclusions 

A new approach to the construction of variational integrators has been developed in this 
work. The new technique is based on the fact that in order to construct the symplectic 
map in the variational integrator, we need only the partial derivatives of the estimation 
of the action integral in a time interval. These derivatives are now functions of the 
integral of the Euler-Lagrange vector, which in the exact case, vanishes. Thus, taking 
an orbit which satisfies the Euler-Lagrange equation in a number of grid points, any 
quadrature used for the calculation of the action integral vanishes and the process is 
dramatically simplified. Experimental tests show that these methods are efficiently 
integrate stiff systems (like the 2-body problem with eccentricity up to 0.99) conserving 
all the benefits of the classical variational integrators. 
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